Astronomy & Astrophysics manuscript no. KHrot 


©ESO 2012 


February 10, 2012 





Impact of orbital motion on the structure and stability of adiabatic 

shocks in colliding wind binaries 

A. Lamberts^ G. Dubus^ G. Lesur\ and S. Fromang^ 

^ UJF-Grenoble 1 / CNRS-INSU, Institut de Planetologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, Grenoble, F-38041, 
France 

2 Laboratoire AIM, CEA/DSM - CNRS - Universite Paris 7, Irfu/Service d'Astrophysique, CEA-Saclay 91191 Gif-sur-Yvette, 
France 

February 10, 2012 

ABSTRACT 

Context. The collision of winds from massive stars in binaries results in the formation of a double-shock structure with observed 
signatures from radio to X-rays. 

Aims. We study the structure and stability of the colliding wind region as it turns into a spiral due to orbital motion. We focus on 
adiabatic winds, where mixing between the two winds is expected to be restricted to the Kelvin-Helmholtz instability. Mixing of the 
Wolf-Rayet wind with hydrogen-rich material is important for dust formation in pin wheel nebulae such as WR 104, where the spiral 
structure has been resolved in infrared. 

Methods. We use the hydrodynamical code RAMSES to solve the equations of hydrodynamics on an adaptive grid. A wide range 
of binary systems with different wind velocities and mass loss rates are studied with 2D simulations. A specific 3D simulation is 
performed to model WR 104. 

Results. Orbital motion leads to the formation of two distinct spiral arms where the Kelvin-Helmholtz instability develops diff'erently. 
We find that the spiral structure is destroyed when there is a large velocity gradient between the winds, unless the collimated wind 
is much faster. We argue that the Kelvin-Helmholtz instability plays a major role in maintaining or not the structure. We discuss the 
consequences for various colliding wind binaries. When stable, there is no straightforward relationship between the spatial step of the 
spiral, the wind velocities, and the orbital period. Our 3D simulation of WR 104 indicates that the colder, well-mixed trailing arm has 
more favourable conditions for dust formation than the leading arm. The single-arm infrared spiral follows more closely the mixing 
map than the density map, suggesting the dust-to-gas ratio may vary between the leading and trailing density spirals. However, the 
density is much lower than what dust formation models require. Including radiative cooling would lead to higher densities, and also 
to thin shell instabilities whose impact on the large structure remains unknown. 
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1. Introduction 

Massive stars possess highly supersonic winds due to radia- 
tion pressure on atomic lines. Wind mass-loss rates range from 
M ^ 10"^Mo yr"^ for O o r B type stars to M ^ IQ-'^Mq yr"^ 
for Wolf-Rayet stars (WR) ( [Puis et al. 12008] ). Most massive stars 
lie in binary systems, where the interaction of the two super- 
sonic stellar winds creates two strong shocks separated by a 
contact discontinuity. The geometry of the colliding wind region 
depends on the momentum flux ratio of the winds ( ,Lebedev &| 
|Myasnikov|199Ql ) 



T] = 



M2V2 
Mivi 



(1) 



where v is the wind velocity. The subscript 1 usually stands for 
the stronger wind, the subscript 2 for the weaker one so that 
T] < I. There are several important observational signatures of 
the colliding wind region. The shoc k-heated gas generates ob- 
servable thermal X-ray e mission (e.g.|Chere pashchuk 1976[|Luol 
eFaL]p^9Ql [Usov|[T9$2l [Stevens et all T992). The presence of 
intra-binary structures caus es variations of the emission line pro- 
files with orbital phase (e.g. [Shore & Brown|1988[|Wiggs & Gi^ 
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1993). The high densities reached in the colliding wind region 
are thought to enable the formation of dust, explaining the large 
infrared emission from binary systems with WR stars ( Williams 
et al.|1987{|Usov|19^T] ), and the formation of spiral structures ex- 
tending to distances up to 300 tim es the binary separation ("pin- 
wheel nebulae", Tuthill et al.|199 9). In some systems, diflTusive 
shock acceleration of particles leads to non-thermal radio emis- 
sion ( |De Becker||2007| ). The radio emission has been resolved 
by long baseline interferometry and s hown to have a morphol - 
ogy changing with orbital phase (e.g. 'Doug herty et aL]|2005 ). 
A new exotic class of colliding wind binaries is gamma-ray bi- 
naries, where the non-thermal emission is thought to arise from 
the interaction of a pulsar relativistic wind with the wind of its 
massive stellar companion ( [Dubus 2006). Interpreting all this 
observational data requires increasingly detailed knowledge of 
the physics of colliding winds, hence numerical simulations, no- 
tably of the large scale regions that can be resolved in radio or 
infrared. 

On large scales, orbital motion is expected to turn the shock 
structure into a spiral, although we will show that this is not 
always true. Orbital motion breaks the symmetry with respect 
to the binary axis and no analytic solution predicts the detailed 
structure of the colliding wind region. Material in the spiral 
is generally thought to behave ballistically, so that the step of 
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the spiral is the wind velocity v tim es the orbital period Porb- 
The wind velocity to use is unclear. 'Tuthil l et al.| ( 2008 ) took 
the speed of the dominant wind v\ (dominant in the sense that 
Mivi > M2V2), whereas Park in & Pittard] (|2008 ) assume that it 
is the slower wind that determines the step of the spiral but fo- 
cus their study on binaries with equal wind velocities. Simple 
dynamical models of the shocked layer have been developed for 
use with radiative transfer codes, assuming that the double shock 
structure is infinit ely thin (thin shell hypothesis) and that the ma - 
terial is balHstic ( [Harries et aLl[2Q04l [Parkin & Pitt^|2QQ8l ). 
The spiral structure is then reproduced at small computational 
cost but this neglects the impact of the pressure, wh ich creates 
a distinction between both arms of the sp iral (Lemaster et al. 
|2Q07l [Pittard|20Q9| [van Marie et"aLl[2QTT] ), the influence of the 
reconfinement of the weaker wind for small 77 (! Lamberts et al.| 
201 1[ ) and the large-scale evolution of instabilities in the collid- 



ing wind region (see below). Up to now no 2D or 3D hydrody- 
namical simulation has modelled a complete step of the spiral. 
We achieve this by using the hydrodynamical code RAMSES 
( jTeyssierj2002| with Adaptive Mesh Refinement (AMR). AMR 
allows large scale simulations to be performed while keeping a 
high enough resolution close to the binary in order to form the 
shocks properly (§2). 



Small sc ale simulations wi thout orbital motion (se e [Stevens] 
let al.|[T992l and references in [Lamberts et al.[[2QTT] hereafter 
Paper I) have shown that several instabilities are at work in col- 
liding wind binaries. Thin shell instabilities occur when cool- 
ing is important so that the shocked zone narrows to a thin 
layer, which is easily perturbed ( [Vishniac|1994| ). They provoke 
strong distortions of the whole colliding region ( Pittard||2QQ9[ 
[van Marie et al.|2QTT ). However, these instabilities are unlikely 
to be dominant in wide binary systems where cooling is in- 
efliicient, th e shocks a diabatic, and the colliding wind region 
wider ( Stevens et al.[|199 2). In this case, the velocity diff'er- 
ence between both winds triggers the Kelvin-Helmholtz insta- 
bility (KHI) at the contact discontinuity (Paper I). Incl uding 
orbital motion has led to contradictory results. Lemaster et al. 
(2007 ) found that eddies develop even when the winds are com- 
pletely identical, because orbital motion introduces a velocity 
diff'erence. [Pitta rd ( 2009 ) found no eddies in a simulation with 
a similar setup, van Marie et al.| (2011 ) also found no eddies, 
although their simulation has an initial non zero velocity diff'er- 
ence p = Vi/v2 = 3/4, and argued that orbital motion stabilises 
the KHI. Larger velocity diff'erences between the winds have not 
been investigated. We performed a set of 2D simulations to study 
how orbital motion changes the shock structure and the develop- 
ment of the KHI, focusing exclusively on binaries with adiabatic 
shocks (§3). The size of the spiral step and the stability of the 
spiral on large scales are addressed in §4. 

WR 104 is an example of a long orbital period binary system 
where the collision between the winds of the WR and its early- 
type companion is expected to be close to adiabatic (see §5). The 
infrared emission is very well matched by an Archimedean spi- 
ral although its brightest point is shifted by 13 milli- arc seconds 
from its centre, possibly because dust formation is inhibited 
closer in dTuthill et al.|2QQ8l hereafter T2008). The WR wind is 
hostile to dust formation due to i ts high temperature, low den- 
sity, and absence of hydrogen ( [Cherchneff" et al. 2000[ ). The 
wind collision region is more favourable, providing high den- 
sities, shielding from the UV radiation of the WR star, and the 
possibility of mixing with hydrogen from the companion star 
( [Marchenko"&M off^at 2007 ). We carried out 2D and 3D hydro- 
dynamical simulations using the parameters of WR 104 to in- 



vestigate these questions (§5). We then relate all our results to 
observations (§6). 

2. Numerical Simulations 

2.1. Equations 

We use the hydrodynamical code RAMSES for our simulations 
(jTeyssier 2002 ). This code uses a second order Godunov method 
to solve the equations of hydrodynamics 



dt 



^ + V • (pv) = 
ot 

-h V • (pvv) -h VP = 



dE 

~dt 



-h V • [\{E -h P)] = 



(2) 
(3) 
(4) 



where p is the density, v the velocity, and P the pressure of the 
gas. The total energy density E is given by 



^ 1 2 P 

2 (r - 1) 

y is the adiabatic index, set to 5/3 to model adiabatic flows. 



(5) 



2.2. Numerical parameters 

We use the MinMod slope limiter together with the exact 
Riemann solver (§3-4) or the HLLC (§5) Riemann solver to 
avoid numerical quenching of instabilities. We perform 2D and 
3D simulations on a Cartesian grid with outflow boundary condi- 
tions. We use AMR which enables to locally increase the spatial 
resolution according to the properties of the flow. We base the 
refinement criterion on velocity gradients. In section §3 we per- 
form small scale simulations where the size of the computational 
domain is l^ox = 40a, with a the binary separation. We have 
Hx = 64 for the resolution of the coarse (unrefined) grid and use 
7 levels of refinement. This gives an equivalent resolution which 
is at le ast two times better than in former s tudies ( [Lemaster et al.[ 
'2OOT; Pittard 2009; van Ma rie etaL|2011| ). In section §4 we per- 
form larger scale simulations where Itox = 400a , the coarse grid 
also has Ux = 64 but we use up to 9 levels of refinement. In some 
cases we adapt the size of our grid to larger or smaller values to 
model a complete step of the spiral. 

2.3. Generation of the winds 

To simulate the winds, we keep the same method as used in 
Paper I, which was largely inspired by [Lemaster et al.[p007 ). 
Around each star, we create a wind by imposing a given density, 
pressure, and velocity profile in a spherical zone called mask. 
The masks are reset to their initial values at each time step to 
create steady winds. We add two passive scalars ^1 and ^2 to 
distinguish both winds and to quantify mixing. We initialise the 
passive scalars in the masks; their evolution is determined by 

dt 



+ V-(p^/v) = /=1,2 



(6) 



In the free wind of the first star = 1 and ^2 = 0, in the second 
wind it is the other way round. In the shocked zone both scalars 
have an intermediate value which accounts for the mixing of the 
winds. The rotation of the stars is clockwise in the figures, their 
positions are updated using a leapfrog method. For each simu- 
lation, the input parameters are the mass M, mass loss rate M, 
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wind velocity v (which we suppose to be constant), and Mach 
number At at r = a for each star. The exact value of the Mach 
number does not matter for the colliding wind region, as long as 
it is high enough that pressure terms can be neglected (Paper I), 
which is the case for massive star winds. Here, the Mach num- 
bers of both winds are set to 30. In all our simulations, the star 
with the highest momentum flux is considered as the first star. 
We will refer to its wind as the stronger wind. The values of the 
parameters of the winds in the simulations are given in the table 
in Appendix B. Both stars have a mass of 15 M© and the binary 
separation a is 1 AU. The corresponding orbital period Porb is 
0.18 yr (67 days). The orbital velocity of the stars is Vorb = 81 
km s"^ We only study circular orbits. 

2.4. 2D and 3D simulations 

We perform our 2D simulations in the orbital plane of the bi- 
nary. We thus model the cylindrical (r, 6) plane instead of the 
(r, z) plane as classically done (e.g. [Stevens et al.|1992[|Brighenti| 



& D'Ercolep995l [Pittard et al.1|2()06| ). This impHes the density 
evolves oc instead of oc in a spherical geometry. For a 
given J] the structure of the colliding wind binary is thus diff'er- 
ent in 2D and 3D. However, as discussed in Paper I, the mapping 
^f^3D V2D captures most of the 3D structure in the 2D simu- 
lations. This point is re-discussed in §5 where we compare the 
results of a 2D simulation of WR 104 with a full 3D simulation 
including orbital motion. A major advantage of our 2D setup is 
the possibility to implement orbital motion for a modest compu- 
tational cost, enabling the study of the flow structure up to scales 
currently inaccessible to full 3D calculations. 

3. Impact of orbital motion on the shock arms 

We carried out simulations of adiabatic colliding winds identi- 
cal to those carried out in paper I except that they now include 
orbital motion to study its impact on the shock structure and 
development of the KHI. The simulations explore t] = 1 and 
T] = 0.0625 for difl'erent velocity ratios jS = Vi/v2 = 1,2,20, 
all in a box of size 8<2. Briefly, the results without orbital mo- 
tion were that (1) no instability is seen when 13=1; (2) for JS > 2, 
the instabilities afl'ect the position of the contact discontinuity, 
for j3 = 20 the KHI also afl'ects the shocks positions; (3) for 
T] = 0.0625 the instabilities remain confined to the weaker wind. 
We present first the results of the simulations for j3= 1 , where the 
KHI instability may be triggered (or not) by orbital shear (§3.1). 
We then discuss the simulations with /3 1. In these cases, the 
dominant wind is slower and much denser than the weaker wind. 
For T] = I, there is no diff'erence between simulations where 
P = Bsindp= l/B. 

A view of the overall colliding wind structure is given in 
Fig. [T] We define the leading arm as the arm preceding the sec- 
ond star, with respect to orbital motion (clockwise motion). The 
trailing arm is the second part of the spiral. Note that there is 
no dominant wind when 7/ = 1 so that the definition of lead- 
ing/trailing is degenerate in this case. (The definition also has no 
link with the definition commonly used in galactic dynamics.) 
In each arm there is a shock in the wind from the first star and 
a shock in the wind from the second star, separated by a contact 
discontinuity. In 2D simulations, when //< 0.25 the second wind 
is confined by the intersection of the shocks. In 3D simulations 
this occurs for rj ^ 0.06 (Paper I). 

Close to the binary, relative motion of the stars creates an 
"aberration" [Parkin & Pitt ard ( 2008 ) of the shocked zone |Parkin| 
& Pittard ( [2008 ) introduce the skew angle yu which measures the 




Fig. 1. Geometry of a colliding wind binary including orbital 
motion ({t]=0. 065^=0.5}). The stars are shown by the black cir- 
cles. The spiral structure has a leading and trailing arm. Each 
arm is composed of one shock from each wind and a contact 
discontinuity. In this case both shocks from the wind of the sec- 
ond star intersect, totally confining the second wind. 



skew angle 77=0.5 Vi = 0.05v2 



0.5 



— 0.0 




Fig. 2. Determination of the skew angle fi for {rj = 0.5, /3 = 0.05}. 
The dashed line showes the theoretical position of the contact 
discontinuity for yu = (no orbital motion). The solid line fits 
the contact discontinuity in the simultion. 

off'set between the line of centers of the stars and the symmetry 
axis of the shocked region. It is given by 



Vorb 

tanyu = 

V 



(7) 



where v is taken as the speed of the slowest wind. This angle 
remains small unless the velocities of the winds are strongly re- 
duced or orbital motion becomes important We measured /j in 
our simulations by finding the best fit of the an alytic position of 
the contact discontinuity derived by [Canto eT al. (1996 ). An ex- 
ample is given on Fig.|2] We measured yu ^ 22° in the simulation 
{rj = 0.0625,yS = 0.05} while the theory predicts fi = 21°. The 
simulation of WR 104 (see §5) gives yu ^ 9° while the theory 
predicts yu = 8°. yu is too small in our other simulations to allow 
correct measurements. 



3.1. Simulations with p=\ 

Fig. [3] shows the density, velocity, and mixing map for a sim- 
ulation with identical winds [q = I, ft = 1}. We determine 
the mixing by the product of the passive scalars ^1 x ^2- The 
free (unshocked) winds correspond to the low density parts at 
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Fig. 3. Small scale simulation: density, velocity, and mixing for a simulation with identical winds (t/ = = 1). The density is given 
in g cm"^, the velocity is in km s"^ the mixing of the winds is a dimensionless variable. The length scale is the binary separation a. 



the top and bottom. The denser parts are the shocked winds. 
The radial inhomogeneities visible in the unshocked wind re- 
gion of the velocity map is a numerical artefact: it corresponds 
to minute anisotropics in the stellar wind due to the finite size of 
the masks. As shown in Appendix A the orbital period is much 
longer than the local shear timescale. The Coriolis force does 
not impact on the development of the KHI. However, the veloc- 
ity map shows that a ^ 20% velocity diff'erence develops in each 
arm at a distance ^ 20a from the binary. Orbital motion makes 
the shocked material leading the contact discontinuity (red re- 
gion on the velocity map) accelerate in the lower density free 
wind region while the shocked wind trailing the contact discon- 
tinuity (in green) moves into the denser, shocked material of the 
other wind. The resulting velocity diff'erence is sufficient to trig- 
ger the KHI even though the original wind speeds are equal. The 
instability is clearly present in the mixing map. We therefore 
confirm the results of |Lemaster et al.| ( [2007| ) that orbital motion 
triggers the KHI even when the winds are identical. We also ob- 
serve, like they do, an artificial enhancement of the instabilities 
when the shocks align with the grid. Our simulations contradict 
the results from v an Marie et al.| ( |20TT] ) who find no KHI. Their 
simulations are performed with a Lax Friedrich Riemann solver. 
We ran a test simulation with the Lax Friedrich Riemann solver 
and observed no development of the KHI either, due to the im- 
portant numerical diff'usivity. 



Fig. |4] shows the density, velocity, and mixing for {ri = 

0. 0625,y5 = 1}. Because of the low value of there is a re- 
confinement shock (Paper I) behind the second star. The various 
discontinuities are indicated in the velocity map, to be compared 
with the simpler geometry shown previously in Fig.[T] Our simu- 
lations without orbital motion showed no KHI because the initial 
velocities are identical. Here, as in the 7/ = 1 case, orbital motion 
leads to velocity shear and mixing at the contact discontinuity. 
The KHI is confined to narrow regions close to the discontinu- 
ity because of the low r] and not because of orbital motion: we 
had found the same behaviour in the models explored in paper 

1. We also see that complex velocity structures arise in the col- 
liding wind region even in this a priori simple case where both 
winds have the same velocity, highlighting the possible difficul- 
ties in interpreting spectral line features arising from this region 
without guidance from numerical simulations. 



3.2. Simulations with JS 1 

Fig. [5] shows the density and mixing for simulations with {7/ = 
= 2,20} and {77 = 0.0625,yS = 0.05,0.5,2,20}. These maps 
show the impact of rotation on arm geometry and the develop- 
ment of the instabilities. 

The leading and trailing arms become markedly diff'erent 
when the velocity diff'erence increases, even when 77=1. The 
shocked zone preceded by the unshocked wind with the higher 
velocity and lower density is larger than the zone preceded by the 
lower velocity, higher density wind. The latter shocked zone is 
compressed by the high velocity wind into a high density region 
fvan Mar ie et al.|2()TT] ). We verified that, as expected if this ex- 
planation is correct, for > 1 the compressed arm is the trailing 
arm while for < 1 the compressed arm is the leading arm (see 
Fig. |5]). For - 20, compression results in a rim of puff'ed up 
matter where the density increases by two orders of magnitude 
with respect to the simulation with /3 = I. The diff'erentiation 
of both arms is independent of the instabilities in the winds but 
plays a role in their development. 

The KHI starts similarly in both arms close to the binary 
major axis as the velocity diff'erence and density jump across the 
contact discontinuity are the same in both arms. The symmetry 
between both arms is broken as the flow moves outwards. The 
compression of the shocked zone in the narrower arm results in a 
thin mixed zone with small scale structures, whereas the eddies 
are stretched out in the wider arm. Mixing covers a larger area 
in the wider zone. This is not just a geometrical efl'ect. We show 
in Appendix A that when media have difl'erent densities (a ^ 0, 
see Eq. A. 10) mixing by the KHI occurs preferentially in the 
least dense medium. Both velocity and density profiles thus play 
a role in the development of the KHI in colliding wind binaries. 
They both impact the large scale outcome of the spiral structure 
as will be shown in ^ 

We checked that the mixing is physical and not numerical in 
the narrow arm. If numerical, the physical size of the mixed zone 
increases with decreasing resolution. If mixing follows from in- 
stabilities, the physical size of the mixed zone is constant with 
resolution. We measure the width of the narrow mixed zone 
across the contact discontinuity at the upper edge of the simu- 
lation box for the case {rj = 0.0625, = 2}. The limit of the zone 
is determined by ^1 x ^2 = 0.01. For the simulations with the 
highest resolution (7 levels of refinement), we measure a width 
of 0.87a while we measure a width of 0.93a in a simulation with 
6 levels of refinement. We conclude that numerical diff'usion has 
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Fig. 4. Small scale simulation: same as Fig.|3]for{77 = 0.0625,yS= 1} 
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a limited impact and that the mixing in our high resolution sim- 
ulations is robust. 



4. Formation of a spiral structure 

We now consider the large-scale evolution of the previous sim- 
ulations in a box of size /^^x = 400a. Density and mixing 
maps are shown for 7/ = 1 in Fig.[6]and for r] - 0.0625 in Fig.[7[ 
with increasing from top to bottom in both figures. The spa- 
tial scale is the same in all plots except for the top two panels of 
Fig.|7] where we reduced the size of the domain to avoid unnec- 
essary computational costs. The diff'erent behaviour of mixing in 
both arms discussed in §3 persists on the larger scales, eventually 
causing both contact discontinuities to merge into one single spi- 
ral in simulations with 77 = 0.0625 (Fig. [7]). We cannot exclude 
this merger results from numerical artefacts due to the use of a 
cartesian grid to describe an inherently spherical phenomenon. 
Merger should still occur naturally from inhomogeneities in the 
winds. Figures [SfT] also show that the colliding wind region does 
not always turn into a stable spiral and that, for given 7/ the ap- 
pearance depends strongly on the velocity ratio fi. We discuss 
below the step size that we measure when a steady spiral forms 
before addressing the issue of the stability of the pattern. 



4. t. The step of the spiral 

We found that, when a stable spiral structure is formed, an 
Archimedean spiral with a step si ze S provides a good fit to 
the results of our simulations. As in |Pittard| (|2009|), we find that 
the fit with an Archimedean spiral is not perfect at the apex. 
However, the deviation is small and limited to a region ^ 10a. 
Mixing follows closely the contact discontinuities in the arms 
so we used the mixing maps to trace the spiral. In fact, the 
spiral is not always clearly apparent in the density maps (e.g. 
{rj = 0.0625,y5 = 1} in Fig.|7]), especially when a complex flow 
is established by the presence of a reconfinement shock behind 
the weaker star (Fig. 1). The fitted S for various values of t] and 
is compared to the theoretical estimate S 1 in Fig. ^ Si as- 
sumes the velocity of the stronger wind controls the structure 
scale so that Si = Porh^i (e.g. T2008). When vi = V2, there 
is no ambiguity in the velocity that sets the step size and we 
verify that, in this case, S = for all t]. This also rules out 
any significant numerical issue with the way the spiral develops. 
There are significant deviations from S i in all the other cases, 
except when 7] <^ 1 i.e. when the first wind largely dominates 
momentum balance. For more balanced ratios 7/, the spiral step 



is smaller than expected when the weaker wind is slower than 
the stronger wind, and vice- versa when the weaker wind is the 
faste st. Using the slowest wind speed instead of vi (e.g. |Parkin| 
|& Pi ttard 2008 ) does not work better. The results do not suggest 
a straightforward analytical correction using 77 and (3 that could 
be used to interpret observations of pinwheel nebulae without 
requiring hydrodynamical simulations. 



4.2. Stability of the spiral structure 

The complete disruption of the spiral structure for large veloc- 
ity ratios (Figs. ]Sff\ was unexpected. The breakdown can be 
traced to the KHI: we found that the structure was stabilised 
when we tested a setup with a Lax-Friedrich Riemann solver and 
a smaller resolution, in which case the KHI is artificially sup- 
pressed (Paper I). The amplitude of the KHI appears to destroy 
the spiral when strong velocity gradients are present, resulting 
in widespread turbulence and important mixing throughout the 
domain. Curiously, for t] = 0.625, the structure is unstable when 
vi = 20v2 while it is stable for the opposite velocity gradient 
vi = 0.05v2: the strong density gradient (M1/M2 = 320) is also 
likely to play a role in the stability. Indeed, we measure a ^ 0.95 
in the mixed region at r ^ 50a, which implies a strong reduction 
of the KHI growth rate. 

Table 1 summarises the presence or absence of spiral struc- 
tures, for all the simulations we performed. We compared these 
results against the KHI growth rate, normalised by the rate at 
which eddies propagate (see Appendix A for details). Although 
the growth rate gives no indication on the saturation in the non- 
linear regime, we suspect that the spiral is destabilized most 
easily when the eddies grow quickly before propagating further 
away. 

Figure [9] shows the normalised growth rate (see Eq. A. 12) 
as a function of j3 for t] = 1,0.5, and 0.0625. For 77 = 1 the 
curve is symmetric. The KHI does not develop when there is 
no velocity diff'erence, peaks at /^max - 2 and drops for higher 
values of j3 because the density gradient dampens the growth 
rate. The symmetry with respect to 13 = 1 is broken when 77 < 1 : 
the normalised growth rate is weaker f or J3 < 1 and stronger for 
j3 > I. The lower the value of 7/, the stronger this asymmetry. 
Hence, stable structures are expected nearyS = 1, for/5 » 1, and 
for 13 ^ 1. In addition, when 7/ 9^ 1, structures with 13 < 1 should 
be more stable than with j3 > I. 

The results of the simulations are in qualitative agreement 
with these expectations. When 7/ = 1, there is a spiral for 13 = 
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Fig. 6. Large scale simulation : density (left column) and mixing 
(right column) for simulations with {t] = l,j3 = 1,2,20} (from 
top to bottom) in the large boxes. The length scale is the binary 
separation a. The mixing of the winds is a dimensionless vari- 
able. 



1,2,4 but not for = 8, 20 (Tab. 1), which is consistent with 
the faster growth of the KHI when increases from 1 . However, 
the transition from stable to unstable spirals occurs further away 
than expected from Fig[9]05max - 2). Also, we were not able to 
recover a stable final structure for very high /3 (or, equivalently in 
the case r] = I, very low /?), up to/5 = 200 (Fig.|9]). Tests at higher 
j3 are computationally too expensive (and may not have much 
astrophysical relevance). For t] = 0.5, we find that the spiral is 
maintained for values close to JS = I but is quickly destroyed for 
higher/lower values of jS. In this case, a stable spiral is recovered 
when /3 < 0.01, consistent with the lower growth rate, while 
the spiral remains destroyed for the symmetric value of = 20 
(higher growth rate). We observe a similar behaviour for r] = 
0.0625. Stabilisation is possible for a higher /3 = 0.05, which is 
consistent with the lower growth rate of the instability for < 1 
as T] decreases. We conclude the presence of a spiral depends on 
T] and in a way that is consistent with having stable structures 
for near-equal velocity winds (vi ^ V2) or when the weaker wind 
is much faster (MiVi > M2V2 and V2 » vi). 



5. The pinwheel nebula WR 104 

WR 104 is a binary composed of an early type star and a WR star. 
The system shows an excess of IR emission related to dust pro- 



g. 5. Small scale simulation: density (left column) and mixing 
right column) for simulations with {7/ = l,yS = 2, 20} (two upper 
rows) and {7/ = 0.0625,yS = 0.05,0.5,2,20} (lower rows). The 
length scale is the binary separation a. 
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Fig. 8. Step of the spiral S/Si as a function of JS for 7] =1 (dia- 
monds), 0.5 (diagonal crosses), and 0.0625 (crosses). The sym- 
bols are respectively joined by a solid, dashed, and dash-dotted 
line for easier identification. 

Table 1. Presence (S) or absence (X) of a spiral structure in sim- 
ulations for various wind momentum (ij) and velocity (J3) ratio. 
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Fig. 7. Large scale simulation : same as Fig. ^ for {tj 
0.0625,yS = 0.05, 0.5, 0.1, 1,2, 20} . 
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Fig. 9. Theoretical 2D growth rate of the KHI in colliding wind 
binaries as a function of the velocity ratio JS = Vi/v2 of the 
winds. The solid, dashed and dash-dotted lines correspond to 
T] = 1, 0.5, 0.0625 respectively. 



duction. The IR emission has been resolved into a spiral structure 
with several steps imaged (T2008). The high temperatures and 
low densities in WR winds are difficult to reconcile with dust for- 
mation, which requires a temperature around 1000 K and a num- 
ber density range between 10^ cm"^ and 10^^ cm"^ dCherchneff" 
|& Tielens||1995| ). An additional constraint for dust formation 
arises from the absence of hydrogen in the WR wind, leading 
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Table 2. System parameters for WR 104 



where A 2 x 



s is the emission rate, 



the 



WR 



OB 



V (km s ^) 
M(Moyr-i) 



1200 (a) 
0.8 X 10-5 - 3 X 10-5 (c) 



2000 (b) 
6x10-8 (d) 



(a) [Howarth & Schmutz 1992| (b) estim ate according to spectral type 
( [Harries et al.|2004|i, (c)[Crowther|1997| ( d) using the mass-loss lumi- 



nosity relation by [Howarth & Prinja|1989| 



to uncommon chemical processes ( [Cherchnelf et al.|2000 ). Dust 
production appears closely-related to binarity and the presence 
of dense colliding wind structures: in eccentric systems, such as 
WR 48 or WR 112, dust production is limited to orbital phases 
close to periastron while it is continuous in systems with circular 
orbits. Systems viewed pole-on show an extended spiral struc- 
ture in infrared. WR 104 is the prototype system of these pin- 
wheel nebulae. Our aim is to determine whether a hydrodynam- 
ical model with adiabatic winds reproduces the observed large- 
scale structure of WR 104, study mixing and identify regions 
where dust production may be possible. Detailed modelling of 
dust formation and growth in colliding wind binaries is beyond 
the scope of this study. 

5.1. Simulation parameters 

Table 2 has the wind parameters for the binary system WR 104. 
The characteristics of the companion to the WR star are not well 
constrained ( jvan der Hucht||20()T] ) and, like T2008, we will re- 
fer to the companion star as the "OB" star. The orbital period 
(241.5+0.5 days), eccentricity e < 0.06, inclination (/ < 16°), 
and angular outflow velocity of the spiral 0.28 mas day"^ in WR 
104 were found by fitting an Archimedean spiral to the IR maps 
(T2008). The orbital separation a is about 2.1-2.8 au for a total 
mass of 20-50 M©. We took e = and a = 2.34 au. Given the 
uncertainties on mass loss rate and velocities rj varies between 
0.0125 = 1/80 and 0.0033 = 1/300. Ass uming a constant ve - 
locity for the OB wind and Rqb = 10 R© ( [Harries et al.|2004) , 
the second shock forms at 2.1 Rqb < r < 5.1/?ob depending on 

The shock position can be influenced by additional physical 
processes. The OB wind is accelerated on distances of ^ 2 - 3 
stellar radii and has not necessarily reached its final velocity at 
the shock, which modifies the eff'ective momentum flux ratio of 
the collision. The shock position moves to 2.27?ob < r < 4.77?ob 
if acceleration is taken into account by using the velocity law 
V = Voo(l -RpB /r). Radiative brakin g of the WR wind by the OB 
radiation field ( [Gayley et al.||1997 1 can also play a role in WR 
104 (T2008). A slower WR wind moves the shock away from the 
OB star (up to 12 Rqb if radiative braking is able to stop the WR 
wind completely, which is only marginally possible in WR 104, 
see T2008). The magnitude of both eff'ects, their compensating 
influence, and the uncertainties in the wind parameters did not 
justify including these processes. We adopted constant velocity 
winds and rj = 0.0033 to ease comparison with T2008. 

Radiative cooling can significantly change the shock struc- 
ture. The ratio x of the cooling timescale tcooi over the dynami- 
cal timescale tesc provides an estimate of its importance ( Stevens 
|etaLri992| ) 



number density of the unshocked wind, ksTs = (3/l6)jdmpV^ 
the shock temperature and Cg the associated sound speed. The 
system is adiabatic if ;^ > 1 and isothermal if ;^ «: 1 . We find 
3 < XOB < 15 for the OB star and 0.3 < ;^wr < 1.2 for the WR 
star. The system is at the transition between the two regimes. 
The escape timescale is assumed to be ^ a/Cs in Eq.[8]but could 
be as short as 2.7 - S.IRqb/Cs (increasing;^ by a factor 10-20) 
if one takes the distance fr om the OB star (~ shock curvature 
radius, Stevens et al.|[1992 ). In the following, we neglected ra- 
diative cooling in the energy equation and assumed an adiabatic 
shock. 

The low value of // is challenging for numerical simulations 
(see discussion in paper I). The mask of the star needs to be as 
small as possible so that the shocks can form properly. A mini- 
mum length of 8 computational cells per direction is needed to 
obtain spherical symmetry of the winds. Numerical resolution 
on scales much smaller than a stellar radius (0.05 au) is thus re- 
quired close to the binary. Further away, we need to maintain a 
high resolution in order to properly study the instabilities, while 
following a spiral step requires a box size > 200 au. We carried 
out two complementary simulations: a 3D simulation covering 
scales up to 12a and a 2D simulation to model a whole step of 
the spiral structure. 

We use the large scale 2D simulation to determine the step 
of the spiral and the impact of mixing. As explained in §2.4, 
we use the mapping aJtjsd ^2D to obtain comparable 2D and 
3D results. We took r]2D = 0.0625 to help comparisons with the 
results in §3-4, which is close enough to 772DWR104 - 0.057 de- 
rived from a straight application of the mapping. It is important 
to have the right velocity diff'erence for the Kelvin-Helmholtz in- 
stability. We thus adapted Mwr in order to have t]2d = 0.0625 for 
the 2D simulation. We use a 200a ^ 500 au simulation box with 
Ux = 128 and 12 levels of refinement. This gives an equivalent 
resolution equal to 2^^ ^ 5 x 10^ cells. We use nested grids to 
slowly decrease the maximum allowed resolution away from the 
binary. 

We use the smaller scale 3D simulations for quantitative re- 
sults on the density and temperature in the winds. The 3D sim- 
ulation follows l/8th of an orbit of WR 104 in a 12a ^ 30 au 
simulation box, large enough to see the impact of orbital mo- 
tion. The orbital plane is the mid-plane of the box and the centre 
of mass of the binary is placed in a comer of the box to max- 
imise the use of the simulated volume. We use adaptive mesh 
refinement with a maximal equivalent resolution of 4096^. We 
limit the high resolution to a narrow zone of 3a close to the bi- 
nary where the instabilities develop. It corresponds to the same 
equivalent resolution as in our 2D model. We model only ^ 20 
layers at this high resolution in the z direction and we gradually 
reduce the resolution when going away from the orbital plane. 



5.2. Global structure 



X 



^cool ^B^^s 
^esc 



(8) 



Figure 10 shows the density, velocity, mixing, and temperature in 
the binary orbital plane of the 3D simulation (top row). The bot- 
tom row has the corresponding 2D map on the same scale. The 
comparison confirms the mapping in r] captures adequately the 
3D structure in the 2D simulation. The positions of the shocks 
and contact discontinuity along the line of centres are similar 
in both 2D and 3D simulation and match the analytic solutions 
(Paper I). The opening angle defined by the contact discontinu- 
ities, well traced in the mixing map, is 15 ± 1°. This angle is con- 
sistent with the analytical estimates that have been used ( [Harries 
et al.[2004| T2008). However, the opening angle defined by the 
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Fig. 10. Density (g cm"^), velocity (km s"^), mixing and temperature (K) in the orbital plane of the 3D simulation of WR 104 (top 
row). Corresponding 2D maps on the same scale (bottom row). The length scale is the binary separation a. 
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Fig. 11. Density (g cm ^), velocity (km s ^) and mixing in 2D simulations of WR 104. The length scale is the binary separation. 



location of the shocks is wider in 2D than in 3D, which may 
have some influence on the density structure at larger scales (see 
§5.3). Because of the low 77, there is a reconfinement shock be- 
hind the OB wind at a distance ^ 0.75a in the 3D simulation 
(1.5<2 in the 2D simulation). All of the OB wind is involved in 
the collision and no fraction escapes freely to infinity. 

Material piles up in both arms of the spiral. T2008 suggest 
diff'erent strengths of the shock can change the conditions for 
dust formation in each arm. In the 3D simulation, the Mach num- 
ber of the trailing arm at r ^ 12a is 13% higher than in the lead 
ing arm, in agreement with the results in|L emaster et al.| ( [2007 



The small temperature diff'erence is unlikely to afl'ect dust forma- 
tion. A more significant eff'ect is that compression keeps a hotter 
temperature in the leading arm than in the trailing arm. Material 
in the mixing zone of the trailing arm experiences a tempera- 
ture an order-of-magnitude cooler than in the mixing zone of the 
leading arm. Dust formation may be favoured in this arm, seed- 
ing the spiral structure when the contact discontinuities merge 
farther out (see below). 

The amount of mixing increases with the distance to the bi- 
nary. Integrating in spheres of increasing radii, the ratio of mixed 



material to the total amount of material within r increases from 
^ 0.01% at r = a to 0.4% at r = 10a. These values are constant 
during the last stages of the simulation (lasting ^ O.lPorb), indi- 
cating the development of the instabilities has reached a steady 
state. As can be seen on Fig. 10 and as expected from theory 
(Appendix A), mixing occurs mostly in the lower density re- 
gions of the colliding wind zone. The velocity map shows that 
the velocity is mostly radial and that matter is accelerated on a 
distance of a few times the binary separation. After substraction 
of the radial component, we find the velocity of the flow along 
the spiral in the 3D simulation reaches a maximal value of ^ 800 
km s"^. This corresponds to the low density region in the center 
of the spiral. In the outer regions of the spiral, the velocity along 
the spiral reaches ^ 500 km s"^ 



Figure[TT] shows the 2D simulation on the largest scale (200a 
or about 470 au). A stable spiral structure forms as expected for 
P - 0.6 and 77 = 0.0625 (Tab. 1). The collimated OB wind gener- 
ates a low density spiral bounded on each side by walls of mate- 
rial where the density is -100 times larger. The initially difl'erent 
mixing in both arms blurs at a distance of ^ 50a. The mixing 
zones more or less merge and follow the leading arm, overlap- 
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ping slightly with the density enhancement of the arm. The step 
of the spiral is 1.055'wr where Swr = vwR^orb = 170 au= 77a. 
T2008 assumed S^r to determine a distance of 2.6 kpc from the 
observed step size. The 5% correction to this distance due to the 
intrinsically larger spiral step is smaller than the uncertainty on 
the measured WR velocity and observed angular step size. 

The single-armed spiral observed in infrared is more remi- 
niscent of the mixing region than the double spiral in the density 
map. A double-armed spiral structure, separated by a very under- 
dense region of angular size ^ 27 mas (at 2.6 kpc), would have 
been resolved if the IR emission correlated with density (i.e. for 
a constant gas-to-dust ratio). However, we caution that the width 
of the low density zone may be overestimated in this 2D sim- 
ulation since it is likely to be related to the opening angle of 
the shocks, which Fig. [TO] shows to be wider in 2D than in 3D. 
Including dust radiative transfer is required for a closer compar- 
ison of the observations with the hydrodynamical simulation. 

5.3. Conditions for dust formation 

One criterion for dust formation is a high enough density. 



Cherchneff & Tielens| ( |1995| ) indicate different paths towards 
the formation of amorphous carbon for number densities n 
ranging from 10^ to 10^^ cm"^ and give a detailed study for 
n = lO^^cm"^. This gives p = 1.4 x 10"^^ g cm"^ assuming a 
mean molecular weight fi = 1.4, typical for a ionized WC wind 
( [Stevens et al.||1992| ). Such a density is only present in our 3D 
simulation at the edge of the spiral, up to a distance ^ 2a from 
the WR star. In the 2D simulation, along the walls we find that 
the density drops asp oc r"^ Using this as guidance, we expect 
p oc r"^ in 3D on large scales. The minimum value ^=10^ cm"^ 
considered by [Cherchneff & Tielens (1995]) is reached at r ^ 25 
a at the inner wall of the spiral. This is equivalent to 1/3 of a 
turn along the spiral. The density is too low for dust formation 
beyond this distance so that any dust present far away has been 
advected out. 

Temperature is another criterion. Models show that dust con- 
densation is possible for 1000 K < T < 6000 K and does not vary 
much within that temperature range (Cherchneff et al.,200 0). A 
strong temperature gradient remains between leading and trail- 
ing regions of the shocked zone, with the more compressed ma- 
terial maintaining a very high temperature. Still, the tempera- 
ture falls from ^ 10^ K close to the stars to ^ 10^ K in parts 
of the mixing zone at the outer edge of the 3D simulation box 

10a). The temperature is expected to fall below 6000 K at 
half a turn of the spiral by extrapolating using T oc r~^^^ (from 
P oc p^/^ oc pT and p oc r"^ in the arms). This can be taken as 
an upper limit to the dust condensation distance as this is cal- 
culated in the adiabatic approximation. It is consistent with the 
infrared observations of a quarter-orbit shift between the maxi- 
mum infrared emission and the binary centre. Radiative cooling 
and photoionization heating of the wind would have to be in- 
cluded to have an accurate determination of the temperature and 
of the impact of shielding from the stellar radiation fields. 



6. Discussion 

6.1. Asymmetries due to orbital motion 

Orbital motion breaks the symmetry around the binary axis and 
introduces significant differences from the stationary case with 
adiabatic winds. It causes a velocity difference that triggers the 
KHI even when both winds are strictly identical. It results in dif- 
ferentiation of the two arms flanking the weaker star. The arm 



moving into the densest unshocked wind is compressed, damp- 
ening the KHI while the other arm expands and sees Kelvin- 
Helmholtz eddies of larger size. The density difference between 
the inner cavity and the b racketing wal l s can r each two orders of 
magnitude. According to Parkin et al. (2011 ) radiation pressure 
can have a similar effect, either enhancing or reducing the ini- 
tial difference. [Varricatt et al.| ( |2004| ) modelled the variations in 
emission line profiles of WR 140 by a rotating cone with dense 
edges, allowing them to constrain the opening angle of the col- 
liding wind region. They found a wider opening angle than ex- 
pected using the analytic formula for the opening angle of the 
contact discontinuity ( Canto et al. 1996) with the standard value 
of 7] for this system. Our simulations also show that matter ac- 
cumulates at the shock rather than at the contact discontinuity. 
The observed opening angle thus corresponds to the opening 
angle of the shocks, which is wider than the opening angle of 
the contact discontinuity. This also increases the fraction of the 
WR wind involved in the collision compared to estimates using 
the contact discontinuity. Some of the spectral line features that 
are not explained by models where both arms have equal emis- 
sion (absorption ) could be due to differen ces between leading 
and trailing arm ( [Stevens & Howarth 1999| ). The skew angle that 
we measured in the simulations matches the theoretical value 
given by Eq. [7] However, as we found that the step of the spi- 
ral is mostly determined by the speed of the stronger wind (and 
not the speed of the slower wind), we wonder whether this could 
als o be the case for the ske w angle. This is implicitely assumed 
by [Kenny & Taylor[ ( [2007[ ). Our simulations do not allow us to 
answer this question. 

6.2. To spiral or not to spiral 

The simulations presented here are the first including at least 
one step of the spiral. We have shown that a structure is main- 
tained on these scales when the two winds have nearly equal 
velocities (J3 ^ 1). This is consistent with the observations 
of pinwheel nebulae in several WR -\- O star binaries ([Tuthill 



[eT^ [20061 [Monnier e t al. ''2007; "Mill our et aL[[2009l ), since 
their winds do have comparable velocities. The spiral is desta- 
bilised when the stronger wind has a velocity between 10-50% 
of the weaker wind (Tab. 1). For example, the episodic ejec- 
tion of large amounts of (initially) slow-moving material could 
have temporarily destabilised any spiral structure i n the lumi- 
nous blue variable (LBV) / WR binary HD 5980 ( [Naze et al. 



2007[[Georgiev et al.|201 1\ . Eta Carina may be a case where any 



large-scale structure generated near apastron (when the system 
is closer to being adiabatic) is destroyed because of the destabil- 
ising velocity ratio j3 ^ 1/6 — although this would have to be as- 



sessed agai nst the effects of the high orbital eccentricity ( Parkin 
et al.pOlT ). We expect our results to hold for eccentric orbits if 



the system stays adiabatic. If the system moves from adiabatic to 
radiative along its orbit then thin shell instabilities develop with 
unknown consequences on the large scale structure. 

The spiral is stabilised again when the velocity ratio JS ^ 1. 
Such a situation may occur in gamma-ray binaries composed 
of a young non-accreting pulsar and an early-type star ( [Dubus[ 
|2006 ). In this case, the stellar wind interacts with the ten uous, 
relativistic pulsar wind. [Bosch-Ramon & Barkov ( 2011 ) have 
argued that the KHI would destroy any large-scale structure. 
Assuming our results hold in the relativistic regime, we find that 
a stable spiral can form on large scales if the stellar wind dom- 
inates because /3 ^ I in this situation. The structure is unstable 
if the pulsar wind dominates, pointing to the intriguing possi- 
bility that the interaction may switch from one regime to an- 
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other in gamma-ray binaries with Be companions such as PSR 
B 1259-63. The highly eccentric orbit takes the pulsar close to 
the equatorial disc where the slow-moving stellar outflow dom- 
inates momentum balance ( |Tavani & Arons|1997| ), leading to a 
stable colliding wind region. However, at apastron, the pulsar 
wind may dominate over the radiatively-driven stellar wind and 
be unable to form a stable structure. Strong mixing of the two 
winds leads to rapid Coulomb or bremsstrahlung losses for the 
high energy par ticles, which ha s an impact on the gamma-ray 
emission (Z dziarski et al.|[2QTQ| ). Extended radio emission was 
detected around PSR B 1259-63 near periastron ([Moldon et al. 



201 la). Regular changes in radio morphology with orbital phase 
have been observed in other gamma-ray binaries that are com- 
patible with non-thermal synchrotron emission from a stable col- 



limated pulsar w ind struc ture on scales < VwPo rh ( [Dhawan et al. 
2006| |Rib6 e Tal. 2008; Moldon et aLl|2011b| ). The luminosity 



and frequency of the radiation are probably too low to be able to 
detect the spiral structure on larger scales. 

An example of colliding winds with y5 » 1 also involves pul- 
sars, this time with a low-mass companion ( Roberts |201 1 ). The 
weak stellar wind is overwhelmed by the relativistic wind of the 
recycled millisecond pulsar. No stable spiral is expected in this 
case. Another possible case is eruptive symbiotics like AG Peg, 
HM Sge or V1016 Cyg. These systems are composed of a red 
giant, with a very slow wind 20 km s"^) and a hot compan- 
ion, a white dwarf in nova outburst at the origin of a fast outflow 
of several 1000 km s"^ ( |Girard & WiUson| 19871 ). We expect the 
spiral structure to be destroyed if the hot companion dominates. 
The radio maps of AG Peg have been interpreted assuming a 
stable spiral structure and a reversal in r] with time (Ke nny &| 
Taylor|[2 007 ). In many cases there is probably no time to form 
a spiral, because of the long orbital period compared to the out- 
burst timescale. The radio maps of HM Sge (possible 90 year 
orbit) show a more fragmented emission region tha n expected 
from colliding wind models (Kenny & Taylor 2005 ), possibly 
because of thin shell instabilities triggered by radiative losses. 

Finally, we note that the KHI also intervenes in the bow 
shock structure created by the adiabatic interaction of the wind 
of a fast-moving star with the interstellar medium (e.g. Mira). 
The controlling parameter is the ratio of wind speed to star veloc- 
ity. Much like with spirals, fast growth of the KHI may strongly 
disturb the cometary structure at large distances, as seen in some 
hydrodynamical simulations ( (Wareing et al. 12006) |2007| ). 



6.3. Dust formation in pinwheel nebulae 

Dust formation in WR 104 and other pinwheel nebulae should 
be helped by the mixing with hydrogen-rich material from the 
early-type star that w e observe in the 2D and 3D simulations. 
Williams et al.| ( |2009 j) argued that stronger dust emission in the 
trailing arm would explain better the IR high-resolution images 
of WR 140, and attributed this to density variations. The winds 
have nearly identical velocities but the WR has a mass-loss rate 
10 times higher than its O companion. The O wind is therefore 
bracketed by two high-density regions. Our simulations do not 
suggest very difl'erent densities. However, larger amplitude mix- 
ing is expected in the trailing arm because it propagates in the 
more tenuous O wind, possibly enhancing dust formation in this 
arm. The lower temperature in the trailing arm also helps. Hence 
a difl'erent dust-to-gas ratio between both arms could be an al- 
ternative explanation. High density eddies triggered by the KHI 
in the arms could also be responsible for the observation of IR 
obscuration event s by dust clouds in other WR-hO star systems 
dVeen et al.|1998l ). 



The ofl'set of the peak IR emission in WR 104 is consistent 
with the distance at which we estimated the temperature to fall 
below the dust sublimation temperature. However, the densities 
in the colliding wind region are on the low side compared to what 
dust formation models require. Adiabatic shocks only enhance 
the density by a factor 4 so radiative cooling is required. Close 
to the binary, the WR wind is likely to present some cooling, 
resulting in a thinner and denser shocked layer. |Pittard ( 2009 ) 
show the post-shock density is about 100 times higher in their 
model cwbl, which has strong cooling, than in their adiabatic 
model cwb3. Radiative cooling also decreases the temperature, 
bringing the region where dust condensation is possible closer to 
the binary. Thin shell instabilities can develop when cooling is 
strong, enhancing mixing of the winds. Given the impact of the 
(weaker) KHI in adiabatic colliding winds, thin shell instabilities 
can also be expected to significantly influence the large scale 
structure. Pittard ( 2009 ) have shown that the diff'erenciation of 
the arms remains when thin shell instabilities are present but the 
large scale outcome has not been studied yet. 

Strong cooling is not necessarily present in all WR binaries. 
I Williams et al.| (|2011|) present evidence from long-term IR ob- 
servations of WR 48 for dust production throughout the orbit. 
The stellar winds in this system have similar characteristics than 
in WR 1 04 but the (tentative) orbital period is much longer, 32 
years. [Williams et al.| ( |2011j ) estimate the system to be adia- 
batic with an average x - 1 1 • The value will be even higher 
at apastron in the eccentric orbit (e = 0.6), yet dust formation is 
present. High densities will be much more difficult to reach than 
in WR 104, requiring dust formation at hitherto lower densities 
than have been considered possible. 

7. Conclusion 

We have studied the large scale impact of orbital motion and 
the Kelvin-Helmholtz instability on adiabatic shocks in colliding 
wind binaries. We used hydrodynamical simulations with adap- 
tive mesh refinement to perform the first simulations of complete 
spiral steps. Orbital motion induces diff'erentiation between both 
arms of the spiral. The arm propagating in the higher density 
wind gets compressed while the arm propagating in the lower 
density wind expands. We explain that this is due to a stronger 
growth of the KHI in the wider arm and discuss possible obser- 
vational signatures in spectral lines. We confirm that the KHI 
arises even when both winds have identical speeds. We compute 
the step of the spiral and caution that there can be large diff'er- 
ences with the standard estimates. We discover that the large- 
scale spiral structure is destroyed when the velocity gradient be- 
tween the winds is important. Strong density gradients have a 
stabilizing eff'ect. According to our simulations we can predict 
whether certain types of binaries present an extended spiral or 
not. Systems with stable spirals are those with near-equal ve- 
locity winds and those where the weaker wind is much faster. 
Performing high resolution simulations of pinwheel nebula WR 
104, we show that in an adiabatic model significant mixing of 
the WR wind occurs with the hydrogen-rich wind of the com- 
panion. The temperature drop allows the formation of dust at 
roughly half a step of the spiral, consistent with the spatial off'set 
in peak IR emission. However, the density in those regions falls 
short of the critical density for dust condensation. Including ra- 
diative cooling would lead to higher densities, and also to thin 
shell instabilities. The impact of these instabilities on the diff'er- 
entiation of the two arms and on the spiral structure is unknown: 
resolving the thin shock layer in a large scale simulation remains 
a very challenging numerical problem. 



11 



A. Lamberts et al.: KHI in colliding wind binaries 



Acknowledgements. AL and GD are supported by the European Community 
via contract ERC-StG-200911. Calculations have been performed at CEA on 
the DAPHPC cluster and using HPC resources from GENCI- [CINES] (Grant 
2011046391) 



References 

Bosch-Ramon, V. & Barkov, M. V. 2011, A&A, 535, A20 
Brighenti, R & D'Ercole, A. 1995, MNRAS, 277, 53 
Canto, J., Raga, A. C, & Wilkin, F. R 1996, ApJ, 469, 729 
Cherchneff, L, Le Teuff, Y. H., Williams, R M., & Tielens, A. G. G. M. 2000, 
A&A, 357, 572 

Cherchneff, 1. & Tielens, A. G. G. M. 1995, in lAU Symposium, Vol. 163, Wolf- 
Ray et Stars: Binaries; Colliding Winds; Evolution, ed. K. A. van der Hucht 
& R M. WiUiams, 346 

Cherepashchuk, A. M. 1976, Pis ma Astronomicheskii Zhurnal, 2, 356 

Crowther, R A. 1997, MNRAS, 290, L59 

De Becker, M. 2007, A&A Rev., 14, 171 

Dhawan, V., Mioduszewski, A., & Rupen, M. 2006, in VI Microquasar 

Workshop: Microquasars and Beyond 
Dougherty, S. M., Beasley, A. J., Claussen, M. J., Zauderer, B. A., & 

Bolingbroke, N. J. 2005, ApJ, 623, 447 
Dubus, G. 2006, A&A, 456, 801 

Gayley, K. G., Owocki, S. R, & Cranmer, S. R. 1997, ApJ, 475, 786 
Georgiev, L., Koenigsberger, G., HiUier, D. J., et al. 201 1, AJ, 142, 191 
Girard, T. & Willson, L. A. 1987, A&A, 183, 247 

Harries, T. J., Monnier, J. D., Symington, N. H., & Kurosawa, R. 2004, MNRAS, 
350, 565 

Howarth, I. D. & Prinja, R. K. 1989, ApJS, 69, 527 
Howarth, I. D. & Schmutz, W. 1992, A&A, 261, 503 
Kenny, H. T. & Taylor, A. R. 2005, ApJ, 619, 527 
Kenny, H. T. & Taylor, A. R. 2007, ApJ, 662, 1231 

Lamberts, A., Fromang, S., & Dubus, G. 2011, MNRAS, 418, 2618 (paper I) 

Lebedev, M. G. & Myasnikov, A. V. 1990, Fluid Dynamics, 25, 629 

Lemaster, M. N., Stone, J. M., & Gardiner, T. A. 2007, ApJ, 662, 582 

Luo, D., McCray, R., & Mac Low, M. 1990, ApJ, 362, 267 

Marchenko, S. V. & Moffat, A. F. J. 2007, in Astronomical Society of the Pacific 

Conference Series, Vol. 367, Massive Stars in Interactive Binaries, ed. N. St.- 

Louis & A. F J. Moffat, 213— F 
Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 
Millour, F, Driebe, T., Chesneau, O., et al. 2009, A&A, 506, L49 
Moldon, J., Johnston, S., Ribo, M., Paredes, J. M., & Deller, A. T. 2011a, ApJ, 

732, LIO-F 

Moldon, J., Ribo, M., & Paredes, J. M. 2011b, A&A, 533, L7 
Monnier, J. D., Tuthill, R G., Danchi, W. C, Murphy, N., & Harries, T. J. 2007, 
ApJ, 655, 1033 

Naze, Y., Corcoran, M. R, Koenigsberger, G., & Moffat, A. R J. 2007, ApJ, 658, 
L25 

Parkin, E. R. & Pittard, J. M. 2008, MNRAS, 388, 1047 

Parkin, E. R., Pittard, J. M., Corcoran, M. R, & Hamaguchi, K. 2011, ApJ, 726, 
105 

Pittard, J. M. 2009, MNRAS, 396, 1743 

Pittard, J. M., Dougherty, S. M., Coker, R. F, O'Connor, E., & Bolingbroke, 

N. J. 2006, A&A, 446, 1001 
Puis, J., Vink, J. S., & Najarro, R 2008, A&A Rev., 16, 209 
Ribo, M., Paredes, J. M., Moldon, J., Marti, J., & Massi, M. 2008, A&A, 481, 

17 

Roberts, M. S. E. 2011, in American Institute of Physics Conference Series, 
Vol. 1357, American Institute of Physics Conference Series, ed. M. Burgay, 
N. D'Amico, P. Esposito, A. Pellizzoni, & A. Possenti , 127-130 

Shore, S. N. & Brown, D. N. 1988, ApJ, 334, 1021 

Stevens, I. R., Blondin, J. M., & Pollock, A. M. T. 1992, ApJ, 386, 265 

Stevens, I. R. & Howarth, I. D. 1999, MNRAS, 302, 549 

Tavani, M. & Arons, J. 1997, ApJ, 477, 439 

Teyssier, R. 2002, A&A, 385, 337 

Tuthill, R, Monnier, J., Tanner, A., et al. 2006, Science, 313, 935 
Tuthill, R G., Monnier, J. D., & Danchi, W. C. 1999, Nature, 398, 487 
Tuthill, R G., Monnier, J. D., Lawrance, N., et al. 2008, ApJ, 675, 698 
Usov, V. V. 1991, MNRAS, 252, 49 
Usov, V. V. 1992, ApJ, 389, 635 

van der Hucht, K. A. 2001, New Astronomical Review, 45, 135 
van Marie, A. J., Keppens, R., & Meliani, Z. 2011, A&A, 527, A3-F 
Varricatt, W. R, Williams, R M., & Ashok, N. M. 2004, MNRAS, 351, 1307 
Veen, P. M., van Genderen, A. M., van der Hucht, K. A., et al. 1998, A&A, 329, 
199 

Vishniac, E. T. 1994, ApJ, 428, 186 

Wareing, C. J., Zijlstra, A. A., & O'Brien, T. J. 2007, ApJ, 660, L129 



Wareing, C. J., Zijlstra, A. A., Speck, A. K., et al. 2006, MNRAS, 372, L63 
Wiggs, M. S. & Gies, D. R. 1993, ApJ, 407, 252 

Williams, R M., Marchenko, S. V, Marston, A. R, et al. 2009, MNRAS, 395, 
1749 

Williams, R M., van der Hucht, K. A., & The, R S. 1987, A&A, 182, 91 
Williams, P. M., van der Hucht, K. A., van Wyk, R, et al. 2011, ArXiv e-prints, 
(arXiv:l 11 1.5194) 

Zdziarski, A. A., Neronov, A., & Chernyakova, M. 2010, MNRAS, 403, 1873 



Appendix A: The Kelvin Helmholtz Instability in 
stratified flows 



A.1. Linear theory 



y 
A 



-u 



p+ 
p- 



Fig. A.l. Configuration of the stratified flow 

We work in the incompressible approximation, assuming we 
have a system with a mean profile U = ±U^^. Abov e j = 0, the 
flow has a density and p" for j < (see Fig. A.l ). We neglect 
the Coriolis force since the local shear timescale ts = Ax/AU ~ 
10"^ yr is much sorter than the orbital period tq ~ 10"^ yr. In 
this approximation, the linearised equation of motion reads: 



p— +pU— + VP = 



dt 



dx 







(A.1) 
(A.2) 



In the following, each quantity is Fourier transformed in x and 
t thanks to homogeneity: Q = Qexp[i(ojt - kx)]. Rewriting the 
equation of motions and combining them leads to 



d^Vy - k\ = 0, 

which is solved with two decaying solutions 

= exp(-ky) for y > 0, 
v~ = A~ Qxp(ky) for y < 0, 



(A3) 

(A.4) 
(A.5) 



and A" being two arbitrarily chosen constants which are ad- 
justed by jump conditions at the interface y = \ pressure should 
be continuous and fluid particles should stick to the interface on 
both sides. The pressure condition reads: 



(A.6) 



where cr- = oj+U. The second condition is obtained defining a 
displacement vector ^(x) which follows the interface. By defini- 
tion, a fluid particle located at (x, ^(x) - e) satisfies 



Vy= — =dt^+ Udjc^ = icr^. 



(A.7) 
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Applying this to both side of interface (±6) leads to the jump 
condition 



(A.8) 



Combining ( A.6 ) and ( A.8 ) and looking for non trivial solutions 
gives 

0? + lacokU + {kUf = 0, (A.9) 
where a = (p^ -p")/(p^ +p"). An instability arises whenever 
A = (l-a^)(-k^U^)<0 (A. 10) 

which is always true since -1 < a < 1. The growth rate is 
1/tkhi = V-A = \kU\ Vl - a^. Hence, a density contrast \a\ 
close to 1 strongly dampens the growth rate of the KHI. 

In colliding wind binaries, the density and velocity of both 
winds are related through the momentum flux ratio //. Using 
Eq. [T] and mass conservation for both winds then, far enough 
from the binary so that ri ^ r2, the density ratio is roughly 

P2 

pi 



(A. 11) 



'^adv 

tkhi 



\p-\\^\-a^ _ 7]^l^p\p-\\ 



(A.12) 



A.2. Nonlinear evolution 
A.2.1. 2D Evolution 

In order to investigate the evolution of the KHI in the nonlinear 
regime, we have performed numerical simulations for increasing 
a. The 2D setup is as follows: box size {Ix = '^Jy = 4), resolu- 



tion (1024x256), code: PLUTO (Mignonee t al.|2007| , adiabatic 
equation of state P oc p^/^, background pressure P = I 'm the ini- 
tial state (using units dimensioned to the box length, density, and 
velocity shear). Reflective boundary conditions are enforced in y 
to confine the instability in the simulation box. We always have 
p^ > p" i.e. the densest medium is found where j > 0. 

In addition to that, we follow a passive scalar s, initialised 
with s = 20(y) - 1 where represents the Heaviside function. 
We performed simulations for {a = 0,0.5,0.9,0.99}. Kelvin- 
Helmholtz eddies are clearly present in the density snapshot 
shown in Fig. |A.2| for model a = 0.5. In order to show the dif- 
fusion of the passive scalar as a function of time, we plo t the 
J sdx SiS 3. function of y and t in Fig. 



A.3 



evolution of ^(y, t) 

These results demonstrate that when a 9^ 0, the scalar diffusion 
propagates much less in the denser medium (y > 0) and that dif- 
fusion looks less efficient when \a\ increases, in the sense that the 
region with intermediate values of the scalar s becomes smaller 
when a increases. 



Fig. A.2. Snapshot of the density at t=21 (in dimensionless units) 
for Of = 0.5. 




Fig. A.3. Diffusion of a passive scalar by the KHI. From left to 
right, top to bottom: a = 0, 0.5, 0.9, 0.99. 




Fig. A.4. Diffusion of the passive scalar in the 3D simulations 
with a = (left) and a = 0.9 (right). 



A.3. 3D evolution 

We have performed simulations for a = and 0.9 in 3D to 
compare them to the 2D ones. They are very similar to the 
2D configuration, except for the resolution which is reduced to 
500 X 100 X 100 in order to reduce computational costs. We set 
= = 4.0. 's(y, t) is shown on Fig. |A.4[ The direct comparison 
with the 2D cases indicates that faster diffusion into the more 
tenuous region is still verified in 3D; diffusion is also slightly 
less efficient in 3D. 



Appendix B: Parameters of the simulations 
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Table B.l. Parameters of the simulations 





vi (km s ^) 


V2 (kms ^) 


Ml (lO-'^Mo yr-i) 


M2 (lO-^Mo yr-i) 


spiral? 


SIS, 


{1,1} 


2000 


2000 


1 


1 


S 


1 


{1,2} 


4000 


2000 


0.5 


1 


S 


0.67 


{1,4} 


2000 


500 


0.25 


1 


S 


0.35 


{1,8} 


4000 


500 


0.25 


2 


X 




{1,20} 


40000 


2000 


0.05 


1 


X 




{1,200} 


8000 


40 


0.05 


10 


X 




{0.5,0.01} 


40 


4000 


100 


0.5 


s 


2.5 


{0.5,0.05} 


200 


4000 


40 


1 


X 




{0.5,0.1} 


400 


4000 


20 


1 


X 




{0.5,0.5} 


1000 


2000 


4 


1 


s 


1.3 


{0.5, 1} 


2000 


2000 


2 


1 


s 


1 


{0.5,2} 


4000 


2000 


1 


1 


s 


0.8 


{0.5,8} 


4000 


500 


1 


4 


X 




{0.5,20} 


8000 


400 


.5 


5 


X 




{0.5, 200} 


8000 


40 


.05 


5 


X 




{0.0625,0.05} 


100 


2000 


320 




s 


1.1 


{0.0625,0.1} 


200 


2000 


160 




X 




{0.0625,0.5} 


1000 


2000 


32 




s 


1.04 


{0.0625, 1} 


2000 


2000 


16 




s 


1 


{0.0625,2} 


4000 


2000 


8 




s 


0.9 


{0.0625,4} 


4000 


1000 


4 




s/x 




{0.0625, 8} 


4000 


500 


4 


2 


s 


0.8 


{0.0625,20} 


40000 


2000 


.8 


1 


s 
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